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The Relative Lyapunov Indicators: Theory and 
Application to Dynamical Astronomy 


Zsolt Sandor and Nicolas Maffione 


Abstract A recently introduced chaos detection method, the Relative Lyapunov In¬ 
dicator (RLI) is investigated in the cases of symplectic mappings and continuous 
Hamiltonian systems. It is shown that the RLI is an efficient numerical tool in de¬ 
termining the true nature of individual orbits, and in separating ordered and chaotic 
regions of the phase space of dynamical systems. A comparison between the RLI 
and some other variational indicators are presented, as well as the recent applica¬ 
tions of the RLI to various problems of dynamical astronomy. 


1 Introduction 


One of the most important questions of investigating a dynamical system with i > 1 
degrees of freedom is to identify the ordered or chaotic behaviour of its orbits. If the 
dynamical system is governed by ordered orbits its time evolution is predictable. On 
the contrary, if the dynamical system evolves through chaotic orbits, its long-term 
behaviour cannot be predicted. In this paper we are considering a special class of 
dynamical systems called Hamiltonian systems. The phase space of a Hamiltonian- 
system usually contains both regions for ordered (predictable) and chaotic (unpre¬ 
dictable) motion, therefore the informations about the locations and extent of these 
regions are of high interest in investigating the evolution of such systems. A typical 
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class of Hamiltonian systems are the planetary systems such as the Solar System, or 
extrasolar planetary systems. 

The ordered behaviour of an orbit or trajectoiyj^Jis strongly related to its stability. 
By the term stability we mean that the trajectories are located in a bounded region of 
the phase space. If the region of chaotic motion is not bounded in the phase space, 
the trajectories could leave that domain through chaotic diffusion. In this case the 
chaotic trajectories become unstable. Thus one way to perform stability investiga¬ 
tions of dynamical systems is the detection of ordered or chaotic behaviour of their 
orbits. The stability of the planets or asteroids in the Solar System is an outstanding 
question of dynamical astronomy. The ongoing discovery of exoplanetary systems 
made the stability investigations of planetary systems even more important. 

In recent years there has been a growing interest in development and applica¬ 
tion of different chaos detection methods. Beside the “traditional” tools such as the 
largest Lyapunov Characteristic Exponent (LCE) or Lyapunov Characteristic Num¬ 
ber (LCN; 0) and the frequency analysis ( ll26ll ). several new methods have been 
developed, which can be used to detect the ordered and chaotic nature of individ¬ 
ual orbits, or to separate regions of ordered and chaotic motions in the phase space 
of a dynamical system. These methods are the Fast Lyapunov Indicator (FLI; Q3; 
l20l ), the Orthogonal Fast Lyapunov Indicator (OFLI; fl6l ). the Mean Exponen¬ 
tial Growth factor of Nearby Orbits (MEGNO; Ifl0li9l). the Spectral Distance (SD; 
ll5Tl ), the Smaller ALignment Index (SALI; BTI l43j| ), the Generalized ALingment 
Index (GALI; B4ll45l ). and finally, the Relative Lyapunov Indicator (RLI; 1 3511361 ). 
whose analysis is the main scope of this paper. We note that all of these methods are 
based on the time evolution of the infinitesimally small tangent vector to the orbit, 
which is provided by the variational equations. Thus these chaos detection methods 
can be classified as variational indicators. 

In what follows, after recalling the definition of the RLI, we present its behaviour 
in symplectic mappings and continuous Hamiltonian systems. The efficiency of the 
RLI is presented by a comparative study with the already mentioned variational 
indicators. This paper closes with a chapter presenting the recent applications of the 
RLI. 


1.1 Definition of the RLI 

The ordered or the chaotic nature of a trajectory can be characterized most precisely 
by the calculation of the LCE: 


1 In the case of continuous dynamical systems the trajectory is a continuous curve in the phase 
space given by the points representing the time evolution of an initial state. In discrete dynamical 
systems the set of the discrete points representing the time evolution of the system is called as 
orbit. 
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where x* £ R" is the initial state of the system, or in other words, the starting point 
of the trajectory, and E, (f) is the image of an initial infinitesimally small deviation 
vector g, between two nearby trajectories after time t. The time-evolution of c is 
given by the equations of motions and their linearized equations: 


where Z)f[x(f)] is the Jacobian matrix of the function f: R” —> R" evaluated at x(f) £ 
R", and x : M — > R", : R — > M" are vector-valued functions, too. In Hamiltonian 
systems if /J (x*) = 0, the orbit emanating from the initial state x* is ordered, if 
l) (x*) > 0 it is chaotic. A serious disadvantage of the calculation of the LCE is that 
it can be obtained as a limit, thus its value can only be extrapolated, which makes 
the identification of weakly chaotic orbits unreliable. 

In practice, one calculates only the finite-time approximation of the LCE, called 
the finite-time Lyapunov Indicator (LI): 


M x ,f) = - Jog 


MM 

Moll 


It is obvious that by calculating the LI for short time, the true nature of individual 
orbits cannot be identified. However, the basic features of the phase space of a sys¬ 
tem (the existence and approximate position of regular regions and extended chaotic 
domains) can be discovered very quickly by calculating a large number of Lis for 
short time. Let x be a vector variable taken along a line, which is going through both 
regular and chaotic regions of the phase space of a dynamical system. Then by fixing 
the integration time fi nt , one can calculate the curve L(x.t lnl ). In the case of regular 
regions (KAM tori, islands of stability) the curve L(x.t nn ) varies smoothly, while in 
the case of an extended chaotic region it shows large fluctuations ( ifTTl L However, in 
the case of weak chaos the fluctuations of the curve L(x,fi n t) are n °t large enough to 
decide the true nature of the investigated region. In order to measure the fluctuations 
of the curve of the finite-time LI at x*, we introduce the quantity: 


AL(x*,t) = |L(x* +Ax*,t) — L(x*,f)| , (1) 

which is the difference between the finite-time LI of two neighbouring orbits, and 
Ax* is the distance between the two initial condition vectors. This quantity has been 
introduced and called RLI in ll35l|36 |. Definition <0 contains Ax* as a free parame¬ 
ter, which should be chosen small enough to reflect the local properties of the phase 
space. In our numerical investigations we have experienced that the arbitrary choice 
of Ax* in a quite large interval 11 Ax* 11 £ [10 14 ,10 7 ] does not modify essentially 
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the behaviour of the RLI as a function of the time. For ordered orbits the RLI shows 
linear dependence on | |Ax* 11, while for chaotic orbits the RLI practically is invariant 
with respect to the choice of |Ax* 11. 

Although there is not developed a strict mathematical theory describing the time 
behaviour of the RLI so far, the results of numerical simulations clearly show its 
power in separating ordered and chaotic orbits. An intuitive explanation could be 
that in the case of ordered orbits the time evolution of the two LI curves (L(x* ,t) 
and L(x* + Ax* .t)) practically cannot be distinguished meaning that they converge 
wiht the same (or very similar) rate to the LCE = 0 limit. On the other hand, the 
convergence rate of the LI of two close chaotic orbits (separated in the phase space 
by the vector Ax*) could be very different, which is reflected in the time evolution 
of the RLI. In the next sections of the paper we shall investigate through extensive 
numerical experiments the completely different behaviour of the RLI as a function 
of time in the cases of ordered and chaotic orbits, which makes it a suitable tool of 
chaos detection. 


1.2 Properties of the RLI in chaos detection 

In order to eliminate the high frequency fluctuations of the curve AL(x,t) for a fixed 
x £ R", we suggested the following smoothing 


i ['Mr] 

(AL{x))(t) = - £ AL(x,i-At) , 

T i =1 


where At is the stepsize. In numerical experiments we always use the above 
smoothed value of the RLI. 

The different behaviour of the RLI for ordered and chaotic orbits are first pre¬ 
sented for discrete Hamiltonian systems, such as the following 2D : 

f x\ = X| +x 2 

\ x ' 2 = x 2 — v ■ sin(xi +x 2 ) mod (2 n) , 


and 4D symplectic mapping: 

! x' l = X\ +x 2 

*2 = X 2 - v- sin(xi +x 2 ) - ju • [1 -cos(xi +x 2 +x 3 +x 4 )] 

X3 = X'3 + X4 

X4 = X4 — K ■ sin(x 3 +X4) — n ■ [1 — cos(xi +x 2 +X3 +X4)] mod (2k) , 

where v and K are the non-linearity parameters, and Ll is the coupling parameter of 
the 4D mapping. 

In the case of the 2D symplectic mapping ([2ji the initial conditions of the ordered 
orbit are xi = 2, x 2 = 0, while the initial conditions of the chaotic orbit are xi = 3, 
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Fig. 1 (a) (left): The phase plot of an ordered and a chaotic orbit in the mapping 0: (b) (right): 
The behaviour of the RLI as the function of time for a chaotic orbit (upper curve) and for an ordered 
orbit (lower curve). 
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Fig. 2 (a) (left): The phase plot of an ordered and a chaotic orbit in the mapping 0; (b) (right): 
The behaviour of the RLI as the function of time for a chaotic orbit (upper curve) and for an ordered 
orbit (lower curve). 


X 2 = 0. In both cases v = 0.5. The phase plots of these orbits are shown in Fig. 
|T[a) and the corresponding time behaviour of the RLI is displayed in Fig. |Tj b). In 
the case of the 4D mapping ([0 the following initial conditions are used: a'i = 0.5, 
a 2 = 0, A 3 = 0.5, A ‘4 = 0 for the ordered, and x\ = 3, aa = 0, X 3 = 0.5, X 4 = 0 for 
the chaotic orbit. In both cases the parameters are v = 0.5, K' = 0.1 and jJL = 0.001. 
The projections of these orbits onto the xi — ao plane are shown in Fig. 0a). The 
behaviour of the RLI as the functions of time of the ordered and chaotic orbits are 
plotted in Fig. 0b). Studying Figs. 0b) and 0b) one can see that the RLI for an 
ordered orbit is almost constant. The RLI of a chaotic orbit grows very rapidly, and 
after reaching a maximum value it decreases very slowly. 

The maximum value of the RLI of a chaotic orbit is much higher (in the examples 
shown by 9-10 orders of magnitude) than the almost constant value of the RLI of an 
ordered orbit. It can be seen that by using the RLI, the ordered or chaotic nature of 
orbits can be identified after a few hundred iterations of the investigated mapping. 

A crucial test for a chaos detection method is whether it separates the weakly 
chaotic orbits (also called “sticky” orbits) from the ordered orbits. In order to 
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Fig. 3 (a) (left): The phase plot of a weakly chaotic orbit in the mapping <|4j; (b) (right): The 
behaviour of the RLI as the function of time for the weakly chaotic orbit seen in the left panel 
(upper curve) and for an ordered orbit (lower curve). 
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Fig. 4 (a) (left): The linear dependence of the final value of the RLI on the initial separation for an 
ordered orbit; (b) (right): The final value of the RLI versus the initial separation for a chaotic orbit. 


demonstrate this property of the RLI we used the 4D symplectic mapping and initial 
conditions for an ordered and a weakly chaotic orbit as EL ED: 

{ x' l = x\-\-x ! 2 

x ! 2 = X 2 + (K/2n)sm(2nxi) — ([l / n) sm[2n(xi — xi)] 

■*3 = *3 +*4 

x ' 4 = X 4 + (K/2n)sm(2nxT,) — (/3 /7r) sin[27r(jri — X3)] mod (1) , 

where K is the non-linearity and P is the coupling parameter. According to ED the 
orbit with the parameters K = 3, p =0.1 and with initial conditions xi = 0.55, X 2 = 
0 . 1 , X 3 = 0.62, X 4 = 0.2 is ordered, while the orbit with the same initial conditions 
and parameter K , but with p = 0.3051 is slightly chaotic tending to a very small 
LCE. The projection of the weakly chaotic orbit on the x\ —X 2 plane is shown in Fig. 
[3ja), and it seems to be an ordered orbit on a torus. Using the RLI (Fig. [3jb)) one 
can see that the chaotic nature of this orbit can be detected after about N ~ 5 x 10 6 
iterations. 
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Finally, we should discuss the role of the initial separation between the two orbits, 
which is one of the free parameters of the RLI (the other one is the lenght of the 
time needed to calculate the RLI, as will be mentioned later on). In what follows, 
we give an evidence that 11 Ax* 11 can be chosen arbitrarily from a quite large interval 
[10 14 ,10 7 ]. The smallest value of this interval is due to the finite representation 
of numbers by computers. On the other hand, the largest value of the above interval 
should also be small enough in order to the RLI reflect the local property of the 
phase space around the orbit under study. 

Fig. ([4]) shows the dependence of the RLI (obtained after a fixed number of itera¬ 
tions, which in this particular case was 2 x 10 4 ), on | |Ax* 11 for the 4D ordered orbit 
of mapping <[3]», shown in Fig.[2|a). In Fig.[4]both the horizontal and vertical axes are 
scaled logarithmically. Studying Fig. [4] a) one can see that for the ordered orbit the 
RLI changes linearly with respect to | |Ax* ||. In Fig. |4j b) we display the dependence 
of the RLI (after 2 x 10 4 iterations) on | |Ax* 11 for the chaotic 4D orbit shown in Fig. 
|2ja). One can see that in the chaotic case the final value of the RLI practically does 
not depend on the choice of ||Ax*||. Thus in the above sense the RLI is invariant 
with respect to the choice of the initial separation. This invariance property can be 
explained by the fact that the RLI of a chaotic orbit practically measures the aver¬ 
age width of the oscillation of the LI curve (as a function of time), which does not 
depend heavily on the choice of the initial separation. 


2 A short comparison of the RLI to other methods of chaos 
detection 


In this section we compare the performances of the RLI with the variational indica¬ 
tors mentioned in Sect. [T| the LI, the FLI and the OFLI, the MEGNO, the SD, the 
SALI and the GALI. 

In Sect. |2.1| we analyze the dependence of the RLI on its free parameters and in 
the following one we compare the typical behaviours of the RLI with other tech¬ 
niques in the well-known Henon-Heiles model ETI (hereinafter HH model). In 


Sect. 2.3 we apply the RLI and other indicators to study dynamical systems of 
different complexity: the 4D symplectic mapping <[3]) presented in Sect. |1.2| and a 
rather complex 3D potential resembling a Navarro, Frenk and White triaxial halo 
(hereinafter NFW model; |[49l ). We compare the phase space portraits given by the 
RLI and the other methods to decide whether the results are comparable. In Sect. 
2.4 we briefly discuss the dependence of the RLI on the computing times. 

The Bulirsch-Stoer integrator is used throughout this section. 
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2.1 The dependence of the RLI on the free parameters 

The RLI has two free parameters: a) the initial separation (Ax*) between the basis 
orbit and its “shadow” (Sect. o and b) its threshold (the threshold is a value that 
separates chaotic from regular motion and it is related with the lenght of the time 
needed to calculate the RLI). These free parameters are “user-choice” quantities. 
Thus, it is of interest to study the dependence of the RLI on both of them. 

The following experiments are undertaken on the HH model: 

& = \ {p 2 x+Py) + \ {* 2 +y 2 ) +x 2 -y- *y 3 , 
where Jf? is the Hamiltonian and x,y,p x ,p y are the usual phase space variables. 


2.1.1 The initial separation parameter 

The dependence of the RLI on the initial separation parameter is strongly related 
to the type of orbit under study (see Sect. ED- Therefore, we take four different 
types of orbits with initial conditions located on the line defined by x = p y = 0 and 
y £ [—0.1 : 0.1] and the energy surface E = 0.118, namely, a regular orbit close to 
a stable periodic orbit (r-sp); a quasiperiodic orbit (r-qp); a regular orbit close to 
an unstable periodic orbit (r-up); and a chaotic orbit inside a stochastic layer (c-sl). 
The initial conditions are taken from (9)- The integration time is 10 4 units of time 
(hereinafter u.t.), which is enough time to provide a reliable characterization of the 
orbits and the stepsize of the numerical integration is 0.01. We note that these values 
have been used in the following numerical experiments, too. 

In Fig.[5ja), we present the final values (i.e. the values of the indicator at the end 
of the integration time) of the RLI as a function of the initial separation parameteij^] 
for the orbits introduced earlier. We show that the initial separation parameter does 
not significantly affect the RLI when we apply the indicator to the chaotic orbit “c- 
sl”, but it does when we apply it to the regular orbits “r-sp”, “r-qp” and “r-up” (and 
confirming the results shown in Sect. This dependence of the RLI on its free 
parameter has severe implications in the selection of the threshold. For instance, if 
we start the computation with an initial separation of 10 l4 , the relation shown in 
Fig. |5ja) will indicate that a good candidate for the threshold to distinguish between 
the chaotic orbit “c-sl” (RLI— 0.1) and the regular orbit “r-sp” (RLI~ 10 135 ) can 
be 10- 12 . Then, the orbits with values of the RLI higher than 10 12 will be classified 
as chaotic orbits. However, this choice of the threshold leads to a misclassification 
of the regular orbits “r-up” (RLI— 10 m j and “r-qp” (RLI— 10 l2 ). Furthermore, 
since the correspondence between the RLI and the initial sepa ration parameter for 
the regular orbits of the sample tends to be linear (see Sect. |l.2) , the threshold 10 12 
does not work at all for an initial separation greater than 10~ u . 

2 The values for the parameter have been taken from the interval suggested in Sect. 1.2 
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Fig. 5 (a) (left): The RLI values as a function of the initial separation parameter for the orbits 
“r-sp”, "r-qp”, “r-up” and “c-sl”; (b) (right): The normalized approximation rates for several chaos 
indicators, including the RLI (see text for further details). 


In order to determine a reliable threshold for the RLI, the relationship with the 
initial separation parameter should be done by computing the indicator for a group 
of orbits known to be regular but with some level of instability (e.g. regular orbits 
close to a hyperbolic object such as an unstable periodic orbit, see E7I ). 


2.1.2 The threshold 

In this section we investigate how the RLI and other indicators (listed above) depend 
on their thresholds. For the following experiment in the HH model we have adopted 
a sample of 125751 initial conditions in the region defined by x = 0, y £ [—0.1 : 0.1], 
p y £ [—0.05 : 0.05] and on the energy surface defined by E = 0.118. The thresholds 
of the LI, the RLI, the MEGNO, the SALI, the FLI/OFLI and the GALIs are shown 
in Table[l] where t is the time (see [13)). From here, the initial separation parameter 
will be 10 l2 , unless stated otherwise. The threshold used for the RLI has been 
computed following the remarks discussed in the previous section. 


Table 1 Thresholds for several indicators, including the RLI 


Indicator 

Threshold 

LI 

ln(t)/t 

RLI 

1(T 10 

MEGNO 

2 

SALI 

1(T 4 

FLI/OFLI 

t 

GALL 

t- 1 

GALL 

t- 2 

GALI 4 

t- 4 
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To proceed with the experiment we define the approximation rate as the rate of 
convergence with a final percentage of chaotic orbits. This rate will show a combi¬ 
nation of the reliability of the indicator and the accuracy of the selected threshold 
if the final percentage of chaotic orbits approaches the “true percentage” of chaotic 
orbits in the system. Therefore, as we require a reliable final percentage of chaotic 
orbits, we consider the percentage of chaotic orbits given by the LI by 10 4 u.t.: 
~ 39.92%, the “true percentage” of chaotic orbits in the system. Both the over¬ 
whelming number of papers claiming the reliability of the LI as a chaos indicator 
and the experimental evidence showing that 10 4 u.t. seems to be a reliable conver¬ 
gent time for all the indicators in the experiment (see, for instance. Sect. |2.2| Fig. 
[ 6 ]l supporting this statement. Thus, we normalize the time evolution of the percent¬ 
age of chaotic orbits given by the methods with this “true percentage”. Hence, the 
values of the normalized approximation rates higher than 1 show percentages of 
chaotic orbits higher than the “true percentage”. 

We test the reliability of the thresholds given in Table [T] according to the above 
mentioned rates. The results are indicated in Fig. [5j b). The convergence towards a 
constant rate of the RLI and the FLI/OFLI is faster than that of the other indicators 
of the sample. Despite this rapid tendency towards a constant value for both indi¬ 
cators, the final percentages of chaotic orbits given by the RLI and the FLI/OFLI 
are higher than that of the LI, which means that the values of the rates are above 
1. Nevertheless, this slight difference in the final percentages of chaotic orbits can 
always be fixed with a small adjustment of the corresponding thresholds. Since the 
results for the MEGNO show substantial disagreement between the percentages of 
the chaotic orbits given by the indicator and the LI, a significant empirical adjust¬ 
ment of the MEGNO’s threshold should be made to avoid an overestimation in the 
number of chaotic orbits. The final percentages given by the SALI and the GALIs 
are in perfect agreement with the “true percentage”. However, their tendency to¬ 
wards a stable percentage of chaotic orbits is slower than the one showed by the 
RLI or the FLI/OFLI. 

Thus, among the above mentioned CIs, the RLI and the FLI/OFLI show the best 
approximation rates, i.e. the best combination of the reliability of the indicators and 
the accuracy of their thresholds. For further details on the experiment, refer to na. 

By 10 4 u.t. the threshold taken for the GALI 4 reaches the computer’s precision 
(10 I6 ) and thus, every chaotic orbit lies beyond such precision. Therefore, the last 
yellow point in Fig.[5jb) falls apart from the tendency. 

Now that we have finished studying the importance of a wise selection of the 
free parameters of the RLI, we calibrated the indicator following the suggestions 
mentioned here and continue comparing its performance with other indicators. 


The Relative Lyapunov Indicators: Theory and Application to Dynamical Astronomy 


11 



10 100 300 1000 10000 


10 100 300 600 3000 10000 



Time [u.t.] 


Fig. 6 Behaviours of (a) the MEGNO, (b) the OFLI, (c) the LI, (d) the SALI, (e) the GALI 3 , (f) 
the GALI 4 and (g) the RLI for orbits “r-sp”, “r-qp”, “r-up”, “c-sl” and “c-cs”. The thresholds as 
well as the lines “I” and “II” are included (see text for further details). 


2.2 Expected behaviour of the indicators in the HH model 

In this experiment, done in the HH model, our goal is to compare the typical be¬ 
haviours of the RLI to the other techniques and show its advantages and disadvan¬ 
tages. 

We take the orbits of Sect. |2.1.T| and a chaotic orbit inside the chaotic sea (c-cs) 
with initial conditions on the line defined by x = p y = 0 and y £ [—0.1 : 0.1] and 
the energy surface E = 0.118 (the initial conditions are taken from 0|). The final 
integration time is 10 4 u.t. and the stepsize is 0.01. 

Fig.[6]shows the time evolution curves of several indicators for the five different 
types of orbits introduced at the beginning of the section. Some of the main features 
of a chaos indicator are the speed of convergence and the resolving power. The 
former is the time it takes to distinguish between chaotic and regular motion. In 
order to visualize this quantity, in Fig. [6] we introduce the vertical lines “I” and “II”. 
The first one shows the time after which the orbit “c-cs” is clearly identified as a 
chaotic orbit with the indicator and the second one plays the same role as “I” for the 
orbit “c-sl”. 
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On panel (a) we present the typical behaviours of the MEGNO (see IflOl ). The 
values of the indicator for the three regular orbits tend towards the theoretical 
asymptotic threshold, 2, in different ways (see the right bottom of panel a). The val¬ 
ues for the chaotic orbits increase linearly with time. At ~300 u.t. (see line “I”), the 
separation between orbit “c-cs” and the threshold line is already significant. Hence, 
orbit “c-cs” is clearly identified as a chaotic orbit then. Only 100 u.t. later (line “II”) 
the same happens with orbit “c-sl”. 

On panel (b) we show the time evolution curves of the OFLI ( lfl 6 l l. The values of 
the indicator for two of the regular orbits increase linearly with time (see the thresh¬ 
old and its expression in Table [TJ while the values for the chaotic orbits increase 
exponentially fast with time. This distinction between both tendencies can be made 
at the same times that have been shown by the MEGNO. Besides, orbit “r-sp” has 
an almost constant value because it is very close to a stable periodic orbit. 

On panel (c) we present the LI (see e.g. 1421 1. The values of the indicator for 
the three regular orbits decrease with time (see the threshold and its expresion in 
Table [T} while the values for the chaotic orbits tend towards a constant value which 
depends on the chaoticity of the orbit. The distinction between the regular orbits and 
the chaotic orbit “c-sl” is made by the LI later than by the previous indicators: the 
orbit leaves the linear tendency of the threshold around --600 u.t. (line “II”). 

On panel (d) we present the time evolution curves of the SALI ( BTl ). The values 
of the indicator for the regular orbits “r-qp” and “r-up” oscillate within the interval 
(0,2), while the orbit “r-sp” tends towards 0 following a power law behaviour. The 
chaotic orbits decrease exponentially fast with time. The time needed for the SALI 
to clearly identify the chaotic orbits is the time used by the chaotic orbits to reach 
the threshold (see its value in Tableland lines “I” and “II” in the figure to locate the 
times). The indicator delays making this distinction (in fact, it does so later than the 
LI) because for smaller values of the integration time, the chaotic orbits decrease 
with a power law as the regular orbit “r-sp” does. 

On panels (e) and (f) we show the time evolution curves of the GALI 3 and the 
GALI 4 , respectively (the GALI 2 and the SALI have almost identical behaviours and 
the former is not included). Their theoretical thresholds (Table |T] see l8l |44I[45l[29l 
for further details) yield good estimations of the time needed for the indicators to 
distinguish the chaotic orbits from the regular orbits. The GALI 3 makes this distinc¬ 
tion in the same time as the LI did (see lines “I” and “II”). Once again, the reason 
for this delay is that the GALI 3 decreases with a power law for regular orbits as well 
as for chaotic orbits at the beginning of the integration interval. Nevertheless, the 
higher the order of the GALI, the faster its tendency towards 0 for the chaotic orbits. 
Then, the GALI 4 has registered the best time so far to distinguish the chaotic orbit 
“c-cs”: ~200 u.t. 

On panel (g) we present the time evolution curves of the RLI (Sect. 1.2 1 . The 
values of the indicator for the three regular orbits are in the interva l (10 l:z , 10 10 ) 
according to the initial separation of 10 ~ 12 (see Fig.[5fa), in Sect. 2.1.1 1 . Thus, as 
the value 10 10 (depicted with a dotted line in the figure) that have been selected in 
Sect. |2.1.2| is in the limit of the interval, it is not reliable as a threshold any more. 
Therefore, we selected the value 10 8 for the threshold. The characterization of the 
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regular orbits does not clearly differentiate among them as the MEGNO, the OFLI 
or the SALI. The values for the chaotic orbits increase with time until they reach 
a constant value. On the one hand, orbit “c-cs” is clearly identified as a chaotic 
orbit by ~80 u.t. (line “I”) when the orbit reaches the threshold. This is the fastest 
characterization of the chaotic orbit “c-cs”. On the other hand, the RLI identifies 
orbit “c-sl” as a chaotic orbit around ^850 u.t. (line “II”), which is the slowest 
characterization. 

All the indicators delay in making a reliable characterization of the chaotic orbit 
“c-sl”, which shows that the chaotic orbit “c-cs” has a larger LI than orbit “c-sl”. 

Finally, the characterization of the five representative orbits made by the RLI as 
well as its speed of convergence is similar to the other techniques. Thus, the RLI is 
most welcome to the group of fast variational indicators. 


2.3 Performances of the indicators under different scenarios 

We have seen in the previous section how similar are the performances of the RLI 
and the other fast indicators in the rather simple HH model. Here we will focus on 
experiments in scenarios that are different from the HH model to determine whether 
the RLI is a reliable technique for studying different or more complex systems than 
the HH model. 


2.3.1 The 4D symplectic mapping 


The time evolution curves of the indicators (used in Sect. |2.2| i are not efficient to an¬ 
alyze a large number of orbits. The appropriate way to gather information in these 
cases is in terms of the final values of the methods. Thus, let us now turn our atten¬ 
tion to the study of the resolving power of the techniques using their final values. 
The following study will be conducted in the 4D mapping Q presented in Sect. 


1.2| by adopting different samples of initial conditions and Hr iterations. The ver¬ 
sion of the MEGNO considered here is the MEGNO(2,0), whose threshold value is 
0.5 (see @). 

The large number of iterations used in the experiments deserves a further expla¬ 
nation. In Fig. 7] we present the RLI mappings for 10 3 (left panel), 5 x 10 3 (middle 
panel) and Krfright panel) iterations. The RLI mapping corresponding to 10 3 it¬ 
erations presents a very noisy phase space portrait probably due to a combination 
of a poor election of the initial deviation vectors (see for instance in El) and the 
short number of iterations. It is also clear from the figure that the phase space por¬ 
trait presents a stable picture after 5 x 10 3 iterations. Furthermore, increasing the 
number of iterations helps to resolve very sticky orbits but no further advantage is 
observed. Thus, we iterate the map 10 5 times in order to distinguish the most sticky 
regions. 
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RLI mapping, 10 3 iterations 



-3 -2.5 -2 -1.5 -1 -0.5 0 
x i 


RLI mapping, 5x10 3 iterations 



-3 -2.5 -2 -1.5 -1 -0.5 0 
x i 


RLI mapping, 10 4 iterations 



-3 -2.5 -2 -1.5 -1 -0.5 0 
x i 


Fig. 7 The RLI mapping on gray-scale (color online) plots composed of 10 6 initial conditions, for 
10 3 (left), 5 x 10 3 (middle) and 10 4 (right) iterations. 



Fig. 8 (a) (left): The RLI mapping on gray-scale (color online) plot composed of 10 6 initial con¬ 
ditions, for 10 s iterations (in logarithmic scale); (b) (right): The MEGNO(2,0), the LI and the RLI 
final values for 10 3 initial conditions along the line xo = —3, for 10 5 iterations (the LI and the RLI 
final values are in logarithmic scale). 


Initial conditions inside high-order resonances 

In Fig. | 8 ja) we present the RLI mapping for a region of the phase space corre¬ 
sponding to the 4D mapping and is composed of 10 6 initial conditions. The main 
resonance as well as the high-order resonances are clearly depicted in dark gray 
(i.e. small values of the RLI) while the stochastic layers inside the main stability 
island and the chaotic sea are depicted in light gray (large values of the RLI). We 
show an horizontal line of initial conditions (x\ £ [—7T.0], xj = —3, X 3 = 0.5 and 
X 4 = 0) used in the following experiment to compare the performances of the RLI 
with the mostly used variational indicator, the LI, and with the MEGNO(2,0), which 
is faster than the LI and which is also a reliable indicator. With a diagonal segment 
we depict the initial conditions (xi =xt€ [—1.03,—0.8], X 3 = 0.5 andx 4 = 0) used 
in the experiment of next section. 

In Fig. [ 8 } b), we compare the performances of the RLI, the LI and the MEGNO(2,0) 
on 10 3 equidistant initial conditions lying on the horizontal line that crosses the 
high-order resonances in Fig. [ 8 ja). This figure clearly shows that the RLI unzips 
the hidden structure inside the high-order resonances better than the LI. Further- 
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(a) (b) 




-1 -0.95 -0.9 -0.85 -0.8 


Fig. 9 Sticky region inside the interval (-1.03,-0.8) for 10 5 iterations and for three indicators: (a) 
the LI, (b) the RLI and (c) the MEGNO(2,0). The representative orbits are depicted with points of 
different types. The thresholds are depicted with a dashed line (the LI and the RLI final values are 
in logarithmic scale). See text for details. The figures were taken from m and slightly modified. 


more, the RLI and the MEGNO(2,0) reveal similar structures (the Y-range of the 
MEGNO(2,0) has been centered on the threshold and shortened to amplify the de¬ 
tails of the revealed structure). For further discussions on the experiment, refer to 

urn 

On behalf of the previous experiments the RLI is not only more reliable than the 
LI to reveal small scale structures, but also an accurate indicator to describe a large 
array of initial conditions. 


Sticky orbits 

Sticky orbits are the most difficult type of orbit to characterize by a variational 
indicator. Thus, we further analyze the identification of this type of orbits (Sect, 
to study the performance of the RLI. 

In Fig. [9} we show the sticky region enclosed in the interval (-1.03,-0.8) (and 
depicted earlier in Fig. |8ja) with a diagonal segment) in terms of the final values 
of the same indicators previously used. We also point out the final values of three 
representative orbits, two chaotic orbits (one of them which is sticky chaotic) and a 
regular orbit. The thresholds are also included. 
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In Fig. |9ja), the LI perfectly indentifies the three representative orbits and their 
domains. In Fig. |9j b), some sticky and chaotic orbits share the same RLI final values 
(~ 1(T 3 ' 5 ) which hide the different levels of chaoticity (such is the case of the rep¬ 
resentative chaotic and sticky chaotic orbits, both of them have similar RLI values). 
This is not the case for the other indicator shown in the figure: the MEGNO(2,0) on 
panel (c). The MEGNO(2,0) has completely different final values for the sticky and 
the chaotic orbits. 

The RLI shows a reliable performance revealing the global characteristics of the 
region, such as the regular domain (where we find the representative regular orbit) 
and some small high-order resonances (e.g. x\ = X 2 ~ —0.85). Nevertheless, it does 
not distinguish the sticky from the chaotic orbit in the experiment (see B71 for 
further details). 

In the following section, we follow with the experiments in a scenario of com¬ 
pletely different nature and much more complex than the HH model or the 4D map¬ 
ping. 


2.3.2 The 3D NFW model 

According to the current paradigm of a hierarchically clustering universe, large 
galaxies formed through the accretion and mergers of smaller objects. The imprints 
of such events should be well preserved in the outer stellar haloes where dynamical 
mixing processes are not significant in relatively short times (for instance, many stel¬ 
lar streams have been identified in the outer regions of the Milky Way (22l l'28lPH). 
Furthermore, this galaxy formation paradigm predicts that the centres of the ac¬ 
creted component of stellar haloes should contain the oldest products of accretion 
events 02 in the formation history of the galaxy (such as the Milky Way), and 
therefore, this substructure might be hidden in its inner regions (e.g. close to the So¬ 
lar neighbourhood). However, in order to study the phase space portraits of stars in 
small volumes in the inner regions of the stellar haloes to quantify and classify the 
substructure, we need a model of the Dark Matter (DM) halo that hosts the galaxy. 

In Il33ll34ll the authors introduced a universal density profile for DM haloes (i.e. 
haloes with masses ranging from dwarf galaxy haloes to those of major clusters): 
the NFW profile. However, in Cold Dark Matter (CDM) cosmologies DM haloes 
are not spherical. Furthermore, numerical simulations suggest that their shape vary 
with radius. In |49| the authors have built a triaxial extension of the NFW profile, 
resembling a triaxial DM halo; the corresponding potential is the so-called NFW 
model. Therefore, the NFW model is a triaxial potential used in galactic dynamics 
associated with equilibrium density profiles of DM haloes in CDM cosmologies. 

In a forthcoming paper we study the phase space portraits of stellar particles in¬ 
side solar neighbourhood-like volumes to gain insights about the formation history 
of Milky Way-like galaxies. Hereinafter, we use some results from that investigation 
to demonstrate the reliability of the RLI on a rather complex model. 

The solar neighbourhood volume is a sphere of 2.5 Kpc of radius located at 8 
Kpc from the center of the NFW model (denoted by <J>n in the following equation): 
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<J>N 





where A is the constant: 


^ __ GM200 _ 

In (1 + Cnfw ) — Cnf w / (1 + Cnf w ) 

with G, the gravitational constant, Afoot), the virial mass of the DM halo and 
Cnfw, the concentration parameter used to describe the shape of the density profile. 
Besides, r p follows the relation: 

_ ( r s + r)r e 

r P — , ! 

r s + r e 

where r s is the scale radius dehned by dividing the virial radius of the DM halo by 
Cv/-'tv. The scale radius represents a transition scale between an ellipsoidal and a 
near spherical shape of the <J>n- The r e is an ellipsoidal radius: 

where b/a and c/a are the principal axial ratios with a the major axis and where we 
require a 2 +b 2 +c 2 = 3 (see [[49]). The values of the constants used in the following 
experiment are taken at redshift z = 0 and listed in Table [2] (see fl9l and references 
therein for further details on the model). 


Table 2 Constants used for the <&n potential 


A 4158670.1856267899 
r s 19.044494521343964 
a 1.3258820840000000 
b 0.86264540200000000 
c 0.70560584600000000 


In order to begin the study of the (6 dimensional) regions of interest in the phase 
space of the NFW model, we needed to restrain some of the variables that defined the 
original sample of 22500 initial conditions. We fixed the positions of the particles to 
the centre of the solar neighbourhood. Then, the stellar particles had the following 
positions at the begining of the simulation: x = 8, y = 0 and z = 0 (in [Kpc]). The 
initial velocity in the polar axis (iv) was restrained to the value —250 in [km s -1 ]. 
The energy ( E ) was restrained to the interval (E m b,Eib) with E m i, ~ —195433 the 
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Fig. 10 Gray-scale (color online) plots of the RLI mapping for the velocity surface v z = —250. 
(left): For a fixed integration time of 10 Gyrs. (right): For a fixed number (150) of radial periods. 
The values of the indicator are in logarithmic scale. 


energy of the most bound particle, and En, ~ —59293 (in [M 0 Kpc 2 Gyr 2 ] with 
[M 0 ] the mass in solar mass units) the energy of the least bound particle of the 
sample. The angular momentum ( L z ) was restrained to the interval (—2000,2000) 
in [Kpc km s -1 ]. We integrated the initial conditions for different integration times 
(1 u.t. corresponds to ~1 Gyr). 

In Fig. |T0} we present the phase space portraits given by the RLI for two differ¬ 
ent choices of the integration times. On the left panel, we integrate the orbits for a 
fixed time interval of the order of the Hubble time, i.e. 10 Gyrs. On the right panel, 
we integrate the orbits for a fixed number of radial periods. The radial period of 
the stellar orbits in galactic potentials such as the NFW model scales as ~ £~ 3 / 2 . 
Then, we integrate the orbits for 150 radial periods in order to have a stable portrait 
of the phase space. It is evident that 10 Gyrs (left panel) is not enough to classify 
properly the orbits with the RLI (or any other indicator). Then, on the right panel 
the time interval used was [~ 57,750] Gyrs where 750 Gyrs is enough time to set 
reliable values of the RLI for the least bound particle of the sample. However, the 
most sticky regions are not clearly depicted yet. Therefore, in the following exper¬ 
iment we choose to scale the integration time linearly with the energy of the orbit. 
The linear relation between the computed integration times and the energy overes¬ 
timates the former for the most bound particles. Indeed, the time interval used for 
the experiment was [~ 204,750] Gyrs where the integration times are clearly larger 
than those applied with the ~ £~ 3 / 2 scale. The larger integration times given by 
the linear scale improve the identification of the most sticky regions which helps to 
evaluate the performance of the indicators. 
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On the left panels of Figs. [TT] and [12] we present the gray-scale plots of the final 
values of four different chaos indicators in the (E, 


norm-, L z ) plan^j the RLI and the 


MEGNO; panels (a) and (bl) of Fig. 11 respectively; the OFLI and the I/GALI 3 : 


panels (cl) and (dl) of Fig. 12 respectively. The GALI 3 is inverted in the color- 
scale plots to make the comparison of the portraits with the other indicators easier. 

In many situations it is useful to define a saturation value by which the chaos 
indicator “saturates”. For instance, the OFLI and the GALI 3 for chaotic motion 
increases or decreases exponentially fast, respectively. Then, if the chaotic nature of 
an orbit is well characterized by the OFLI when the indicator reaches 10 16 or by the 
GALI 3 when it reaches the computer’s precision (10~ 16 ), the computation should be 
stopped. Hence, the values 10 16 and 10 16 can be used as saturation values for the 
OFLI and the GALI 3 , respectively. Another example is the MEGNO: the MEGNO 
has an asymptotic value for regular orbits, 2 (see Table]!}, and increases linearly for 
chaotic orbits. Then, if the MEGNO reaches the value 30, the orbit is undoubtedly 
chaotic and it is worthless to continue the computation of the indicator. Then, the 
value 30 can be used as a saturation value for the MEGNO. The time of saturation, 
that is the time by which the indicator saturates, it is a quantity useful in recovering 
the chaoticity levels in the chaotic domain: the smaller the value of the time of 
saturation, the more chaotic the orbit (see SHEzl). Finally, if the indicator saturates 
(i.e. the indicator reaches its saturation value), the integration times will be the times 
of saturation, but, if the indicator does not saturate, the integration times will be the 
final integration times. 


On the right panels of Figs. 11 and 12 we present gray-scale plots of the inte¬ 
gration times used for three of the four indicators above mentioned. On panel (b2) 


Tl| we present the integration times for the MEGNO or MFGNO SY „ and in 
panels (c2) and (d2), the integration times for the OFLI and the GALI 3 , or 
u[ li m/ and GALL,'", respectively. 

The discussion below is not intended to analyze the dynamics of the system, 
which is the aim of a future work. Our goal here is to demonstrate that the perfor¬ 
mance of the RLI in this complex system is as good as the performances of the other 
three wide-spread indicators. 


On the left panels of Figs. 11 and 12 we can clearly see the great level of coinci¬ 
dence among the phase space portraits of the four indicators. The regular component 
is composed of symmetrical structures around the L z axis and the four indicators 
represent these structures with very similar shapes, sizes and shades of gray. This 
shows that the indicators do not only agree in the location, extension and shape of 
the domains of regular motion, but also in the description of these domains. 

The chaotic domain is described equivalently by the four indicators and their 


corresponding times of saturation. For instance, the RLI (Fig. 11 a)) shows that the 
chaoticity (light gray) is inversely proportional to E norm and does not depend on the 
L z . That is, more bound orbits are more chaotic orbits. We arrive at the same (triv¬ 
ial) conclusion with the information given by the integration times of the other three 
indicators (Figs. ]TTjb2),|~i~2}c2) and 1 1 2 [d 2 )). The MEGNO (Fig. 11 (bl)) shows an 


3 The E norm is the normalized energy: (E - E mb ) / {E tb - E mb ). 
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Fig. 11 Gray-scale (color online) plots for the velocity surface v z = —250. (a) The RLI mapping 
(the values of the indicator are in logarithmic scale), (bl) the MEGNO mapping and (b2) the 
MEGNO ra i mapping. 


uniform and almost white color for the chaotic domain. This implies that the satu¬ 
ration value (30) has been reached by these chaotic orbits and thus, no further struc¬ 
ture is revealed. However, the MEGNO va/ (Fig. 1 1 1 |b2)) shows such structure (sur¬ 
rounded by orbits that did not saturate^] This structure revealed by the MEGNO sa r 
shows that the times of saturation are directly proportional to E, wrm or, once again, 
that chaoticity is inversely proportional to E norm and also does not depend on the 
L z . Similar conclusions can be drawn from Fig.[l2]for the OFLI (panel cl) and the 
OFLI jaf (panel c2) and the GALI 3 (panel dl) and GALI"" (panel d2). However, the 
region composed of orbits that have reached the associated saturation values ( 10 16 

4 Remember that the integration time varies with the E„ orm , which explains the transition from dark 
to light gray in the background of the plots on the left side of Figs.[TT]and[T2] Also the time of 
integration is fixed to 0 where there are not initial conditions. 
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Fig. 12 Gray-scale (color online) plots for the velocity surface v z = —250. (cl) The OFLI mapping 
and (c2) the OFLI jaf mapping, (dl) the GALI 3 mapping and (d2) the GALI™' mapping (the final 
values of the indicators are in logarithmic scale). 


and 10 16 for the OFLI and the GALI 3 , respectively) within the interval of integra¬ 
tion is now much extended. On the one hand, this region in Figs. [I2|'cl) and (dl) 
is depicted in an uniform and almost white color and thus, the structure cannot be 
revealed. On the other hand, the times of saturation in Figs. [T2jc2) and (d2) fulfill 
the missing information. 

In the next experiment, we follow two orbits in the NFW model for a time-span 
of 1000 Gyrs (i.e. ~ 77 Hubble times) in order to have convergent final values of 
all the indicators in the study. We compute the time evolution curves of several 
indicators, including the RLI, and present the results for a chaotic and a regular 
orbit (“cha” and “reg”, respectively) in Fig. 13 In order to distinguish efficiently 


the chaotic orbit from the regular one, we proceed as in Sect. |2.2| and use the same 
thresholds used there for the MEGNO, the OFLI, the LI, the SALI and the RLI. 


OFLI, 
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CIS vs RLI {NFW model) 



Time [u.t.] 


Fig. 13 Behaviours of (a) the MEGNO, (b) the OFLI, (c) the LI, (d) the SALI, (e) the GALI 3 , (f) 
the GALI 5 and (g) the RLI for orbits “cha” and “reg”. The thresholds as well as the line “I” are 
included (see text for further details). 


The threshold for the GALI 3 in the NFW model will be the same constant used for 
the SALI (see ( 8 )) because the model is a 3 degree of freedom (hereinafter d.o.f.) 
system. The threshold for the GALI 5 will be t 4 , with t the time (see El for further 
details). In Fig. 13 we mark with the vertical line “I” the time after which the orbit 
“cha” is clearly identified as a chaotic orbit. 

Fig.[l3]shows that the accurate identification of the orbit “cha” as a chaotic orbit 
by the MEGNO (panel a), the OFLI (panel b), the GALI 5 (panel f) and the RLI 
(panel g) is made within a Hubble time (~13 Gyrs). The above mentioned indicators 
show that the orbit will behave as a chaotic orbit within a physical meaningful time- 
span (i.e. the age of the Universe) which is important to understand the dynamics of 
a real system like a galaxy. 

In this section, we present results that support the RLI as a reliable indicator. This 
technique shows a phase space portrait very similar to those shown by the MEGNO, 
the OFLI or the GALI 3 . Furthermore, the RLI identified the chaotic nature of the 
orbit “cha” very fast in the second experiment (the fastest being the OFLI and the 
GALI 5 ), within a Hubble time. These results put the RLI on an equal footing with 
the other fast variational indicators. 
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2.4 A short discussion on the computing times 

The computing times of the indicators are specially crucial for time-consuming 
simulations and their estimation helps for an efficient usage of the computational 
resources. However, such estimation is not an easy task. The computing times de¬ 
pend on a wide variety of factors such as the complexity of the model and that of 
the indicator’s algorithm, its numerical implementation, the hardware, etc. Although 
making a detailed study of the computing times is not our concern (see fl3l for fur¬ 
ther information on the subject), we would like to point out the fact that the easy 
algorithm of the RLI helps for a rather fast computation of the indicator. Further¬ 
more, we are not dealing with the computing times themselves. We are going to 
register the ratios between the computing times of the different techniques and the 
computing time of the Ll^jfor orbits in two of the systems previously studied in the 
chapter, the HH and the NFW models. If the ratio is above 1, then the computing 
time of the indicator is larger than the computing time of the Lis. 

In this experiment, we used the following hardware/software configuration: an 
Intel Core i5 with four cores, CPU at 2.67 GHz, 3 GB of RAM, an OS of 32 bits, 
and the gf ortran compiler of gcc version 4.4.4, without any optimizations. The 
code to compute the indicators is the LP-VIcode, the acronym for “La Plata Vari¬ 
ational Indicators code”. The alpha version of the LP-VIcode was introduced in 
III and currently, it is a fully operational code that computes a suite of ten varia¬ 
tional indicators (see f 6 )^] The initial setup of the LP-VIcode is the following: 
the integration time is 1000 u.t. (or 1000 Gyrs in the NFW model), the step of inte¬ 
gration is 0.001 u.t. (or 1 Myr in the NFW model). 

We registered the ratios for two pairs of orbits and for every indicator. Both pairs 
have one chaotic and one regular orbit. One of the pairs of orbits is located in the HH 
model and the other in the NFW model. The computing time of the Lis for the HH 
model (i.e. 4 Lis) is 0ml5.204s., and for the NFW model (i.e. 6 Lis) is 3m27.703s. 
The energy conservation is AE ~ 10 ! 1 . 

The results are shown in Tabled The value “N” is the number of d.o.f. of the 
system. For instance, in the HH model (second column in Table |3j, we compute 4 
Lis and 3 GALIs: the GALL, the GALI 3 and the GALI 4 . 

The RLI is not the least consuming indicator so far, according to the information 
given in Table [3] The FLI/OFLI, the MEGNO, the SALI and the SD might be more 
desirable options (in that order). Nevertheless, the easy implementation of the RLI 
helps to reduce significantly the computing time of the Lis, and its ratio is not much 
larger than those of the SD or the SALI. 

5 Here, the Lis are the numerical approximations of the spectra of Lyapunov Characteristic Expo¬ 
nents. 

6 Further information on the LP-VIcode can be found at the following url: 

http://www.fcaglp.unlp.edu.ar/LP-VIcode/ 
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Table 3 Ratios of the computing times for several indicators, including the RLI, for the HH and 
the NFW models 


Indicator(s) 

Ratios (HH model) 

Ratios (NFW model) 

(2N-1) GALIs 

-2.026 

-1.024 

(2N) Lis 

1 

1 

RLI 

-0.892 

- 0.443 

SD 

-0.613 

-0.375 

SALI 

-0.582 

-0.372 

MEGNO 

-0.53 

-0.174 

FLI/OFLI 

-0.432 

-0.151 


3 Application of the RLI to planetary systems 

The ongoing discovery of exoplanetary systems is certainly the most rapidly grow¬ 
ing field of astronomy. Up to now the number of known planetary systems is almost 
200. On the other hand, the masses and orbital elements of the planets detected can 
be determined only with uncertainties. Therefore it is of high importance to pro¬ 
vide stability estimates of these systems. Before the launch of the space missions 
devoted to detect terrestrial planets (CoRoT, Kepler), only the massive giant planets 
were discovered mainly by the radial velocity method. One of the major applica¬ 
tions of the RLI was to study the stability of the still hypothetical terrestrial planets 
in planetary systems containing at least one giant planet j37l . Another possible ap¬ 
plication area of the RLI is to study the stability of different orbital solutions of 
resonant exoplanetary systems provided by radial velocity observations. In this re¬ 
view we shortly describe the stability studies done for the resonant planetary system 
HD 73526 l39l . Finally, we present the applicability of the RLI to map the high 
order resonances in the restricted three-body problem, which might have relevance 
when studying the behavior of Kuiper belt objects 03. We note that close relatives 
to our investigations are the works of ED computing the stability maps of the sys¬ 
tem 55 Cancri by using various indicators (LI, SALI and FLI), and of [7] studying 
the dynamical stability of the Kuiper-Belt using the LI indicator. 

In all of the above problems the mean motion resonances (MMR) play an essen¬ 
tial role, therefore we shortly summarize their properties. A MMR occurs between 
two bodies orbiting a more massive body if their orbital periods can approximately 
be expressed as a ratio of two positive integer numbers, 7) /T) = (p + q)/p, where 
7) and 73 denotes the orbital periods of the two bodies, respectively. A MMR can be 
characterized by studying the behaviour of the resonant angle, which in the model 
of the restricted three-body problem for an inner MMR is 

9 = (p + q)X' — pX — qU5, (5) 

while for an outer MMR can be written as 


6 = (p + q)X — pX' — qO. J, 


( 6 ) 
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where A and £U are the mean orbital longitude, and longitude of the pericenter of 
the massless body, while A' is the mean orbital longitude of the massive body. If 
0 librates with a certain amplitude, the two bodies are engulfed in the (p + q) : p 
MMR. In this way almost the same orbital configuration of the bodies involved in 
the given MMR is repeated. Depending on the relative positions of the bodies, this 
configuration can be protective, or can result in unstable orbits, see for more details 
in J32). 


3.1 Stability catalogue of the habitable zones of exoplanetary 
systems 

The main idea behind the stability catalogue was to map the regions of a plane¬ 
tary system that can host dynamically stable terrestrial planets 03 . The dynamical 
stability of a terrestrial planet is one of the strongest requirements for a habitable 
planetary climate. The most important requirement for the habitability of a planet is 
to contain water in liquid phase on its surface. A region around a star, in which an 
Earth-mass planet could be habitable in the above sense is called as the habitable 
zone (HZ), see l24l and ll25l for more details. 


3.1.1 Used models and initial conditions 

The stability of terrestrial planets can be studied by using different approaches: (i) 
by detecting the stable and unstable regions of the parameter space of each exoplan¬ 
etary system separately or (ii) by using stability maps computed in advance for a 
large set of orbital parameters. In the stability catalogue we presented such stabil¬ 
ity maps also showing how to apply them to the exoplanetary systems under study. 
This second approach has the advantage that the stability properties of a terrestrial 
planet can be easily reconsidered when the orbital parameters of the giant planet of 
an exoplanetary system are modified. This is very often the case, since the orbital 
parameters of the giant planets are quite uncertain, and due to the accumulation and 
improvement of the observational data, they are subject to change quite frequently. 
Instead of the re-exploration of the phase space of each individual exoplanetary 
system after possible modification of the orbital parameters of the giant planet, the 
stability properties of the investigated planetary system can be easily re-established 
from the already existing stability maps. These stability maps, which form a stabil¬ 
ity catalogue, can also be used to study the stability properties of the habitable zones 
of known exoplanetary systems. 

The majority of planetary systems which are detected so far consists of a star 
and a giant planet revolving in an eccentric orbit. Therefore, we used a simple dy¬ 
namical model, the elliptic restricted three-body problem in which there are two 
massive bodies (the primaries) moving in elliptic orbits about their common center 
of mass, and a third body of negligible mass moving under their gravitational in- 
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fluence (for details see H46I ). In our particular case the primaries are the star and 
the giant planet, and the third body is a small Earth-like planet, being regarded as 
massless. We note that among the extrasolar planetary systems there is a high rate 
of multiple planet systems. Thus, a more convenient model for the stability maps 
would be the restricted N-body problem (with N-2 giant planets, N > 4). The pres¬ 
ence of additional giant planets certainly enhance the instabilities induced by just 
one massive planet, turning the HZ of the system more unstable. The main source 
of instabilities are the mean motion resonances (MMR) between the massive giant 
planet and the Earth-like planet. Thus, by mapping these resonances, we can find 
the possible regions of the instabilities in the HZs. On the other hand, the dynamical 
model with one giant planet also offers the most convenient way to display the most 
important MMRs as a function of the mass ratio of the star and the giant planet, and 
of the eccentricity of the giant planet. 

In the catalogue of dynamical stability one important quantity is the mass pa¬ 
rameter of the problem fi = mi /(mo + mi), where mo is the stellar, and m\ is 
the planetary mass. The mass parameter has been changed between broad limits 
(10“ 4 — 10 2 ) with various steps of A u, in total the different stability maps have 
been calculated for 23 values of u. The giant planet was placed around the star in 
an elliptic orbit, with semi-major axis a\, eccentricity e\, argument of periastron C0\ 
and mean anomaly M\. The semimajor axis a\ was taken as unit distance a\ = 1 
during all simulations. The eccentricity e\ was changed between 0.0 and 0.5 with 
a stepsize of 5 x 10 3 . The argument of periastron was fixed at C0\ = 0°, while the 
mean anomaly M\ was changed between 0° and 360° with AM\ = 45°. The test 
planet was started in the orbital plane of the giant planet with an initial eccentricity 
e = 0, argument of periastron ft) = 0°, and mean anomaly M = 0°. The semi-major 
axis a of the test planet was changed in two different intervals: (i) for orbits of the 
test planet ‘inside’ the orbit of the giant planet between 0.1 and 0.9 with a step- 
size of Aa = 10 3 and (ii) for ‘outside’ orbits between 1.1 and 4.0 with a stepsize 
of Aa = 3.625 x 10 3 . Further details and the complete catalogue can be found at 
http://astro.elte.hu/exocatalogue/index.html. 


3.1.2 Stability maps 

Due to the very good visibility of the outer MMRs, we first display the case when the 
semi-major axis of the test planet is larger than the giant planet’s semi-major axis 
(a > 1). Additionally, in the stability map shown in Fig. [14] the values /i = 0.001 
and M = 0° were kept fixed. For each value of the RFI a gray shade has been as¬ 
signed. White regions correspond to small RFI values, thus they are very stable. 
The MMRs appear either as light strips in the dark, strongly chaotic regions or as 
the well-known “V”-shaped structures representing the separatrices between reso¬ 
nant and non-resonant motion. The inner regions of the resonances may be lighter 
than the lines of the bounding separatrices indicating regular motion in a protective 
resonance (e.g. the 2:5 MMR). Near the separatrices the motion is always chaotic, 
moreover at some MMRs even the inner part of the resonance is chaotic as indi- 
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Fig. 14 Stability map of the outer MMRs for the Earth-like planets in the elliptic restricted three- 
body problem for ju = 0.001 and M = 0°. White colour denotes ordered motion, light grey strips 
and “V” shapes the different resonances, while black the strongly chaotic regions. 


cated in the case of the 1:3 MMR, for instance. By increasing the giant planet’s 
eccentricity e\ many resonances overlap giving rise to strongly chaotic and thus 
very unstable behaviour. The reason of this phenomenon is that by increasing e\ 
the apocenter distance of the giant planet also increases, therefore the giant planet 
perturbs more strongly the outer test planet. 

In stability maps displayed in Fig. [15] the semi-major axis of the test planet is 
smaller than that of the giant planet (a < 1). The two panels for the mass param¬ 
eter /i = 0.001 show the stability maps for the two different starting positions of 
the giant planet, M\ = 0° (upper panel) and M\ = 180° (lower panel), respectively. 
Between the test planet and the giant planet, a large number of inner MMRs can 
be found, which dominate the stability maps. Inside the resonances the stable or 
chaotic behaviour of the test planet depends on the initial angular positions of the 
two planets. This is clearly visible by comparing the two panels. On the other hand 
the location of the MMRs is not altered, since this depends on the ratio of the semi¬ 
major axes of the two planets. In the lower panel of Fig. 15 several MMRs (5:2, 
5:3, 3:2) are stable, which is not the case in the upper panel of Fig.[l5]This is due 
to the fact that the relative initial positions determine the places of conjunctions 
of the two planets. If they meet regularly near the pericenter of the giant planet, 
the motion of the test planet becomes chaotic, while it can remain regular if the 
conjunctions take place near the apocenter of the giant planet. The effect of the 
initial phase difference between the planets is important, therefore a bunch of sta¬ 
bility maps have been prepared for more initial values of the mean anomaly M\ 
of the giant planet. These stability maps can be found in the online exocatalogue 
(http: //astro . elte . hu/exocatalogue/index. html). By increasing 
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Fig. 15 Stability map of the inner MMRs for the Earth-like planets in the elliptic restricted three- 
body problem for fl = 0.001, M = 0° (upper panel), and M = 180° (lower panel), respectively. 
We note the different character of the 5:2, 7:3, and 3:2 MMRs depending on the orbital positions 
between the test particle, and the perturbing body. White colour denotes ordered motion, light 
grey strips and “V” shapes the different resonances, while black the strongly chaotic regions. We 
displayed the properly scaled HZs of three exoplanetary systems in the stability maps. Studying 
the figures one can conclude that the HZ of the Solar system is stable, the HZs of e Eridani and 
HD 114729 are marginally stable meaning that they are filled with several MMRs. 


the value of the mass parameter Li it becomes clearly visible that the larger mass 
of the giant planet results in stronger perturbations, and therefore more enhanced 
chaotic region. 
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3.1.3 Stability of terrestrial planets in the habitable zones 


In this section, we show how to use the catalogue to determine the stability of hy¬ 
pothetical Earth-like planets in exoplanetary systems. As an example, we consider 
the case of HD10697, where a\ = 2.13 AU, ei = 0.11 and fi = 0.0055. Fig. [16] 
shows a stability map, calculated for /i = 0.005 for inner orbits of the test planet. 
This corresponds to the minimum mass of the giant planet (minimum masses are 
used throughout). The stability of a small planet (starting with e = 0) in the sys¬ 
tem HD10697 can be studied along the line e\ = 0.11. One can see that for small 
semimajor axes, a < 0.33a i = 0.729 AU the parameter space is very stable. When 
a > 0.33a i several resonances appear, among which the most important are the 5:1, 
4:1, 3:1 and 2:1 MMRs. For a > 0.73a| = 1.55 AU, a strongly chaotic region ap¬ 
pears. The classical HZ of this system is between 0.85 and 1.65 AU therefore in Fig. 
16 the scaled classical HZ is located between 0.85/ai = 0.39 and 1.65/ai = 0.77 
(shown as a rectangle, elongated in horizontal direction). One can see that the inner 
part of the classical HZ contains ordered regions, but stripes of certain resonances 
are also present. The outer part of the classical HZ is in the strongly chaotic region. 
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Fig. 16 Stability map for inner orbits of an Earth-like planet, when fd = 0.005 and Mi =0°. The 
line e\ =0.11, corresponds to the system HD 10697. The scaled HZ of this system is between 
0.39 < a < 0.77. Its inner part is stable (containing only a few weakly chaotic MMRs), while the 
outer part is in the strongly chaotic regions. 
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Fig. 17 Stability maps calculated around the orbital elements of ED- The structure of the stabil¬ 
ity maps indicates that the orbital elements (marked by +) are in a weakly chaotic region. Here 
white colour refers to ordered, lighter grey shades to weakly chaotic, and darker shades to unstable 
regions. 


3.2 Stability of resonant exoplanetary systems 

A significant amount of multiple extrasolar planetary systems contain pairs of giant 
planets which are orbiting in MMRs. The in situ formation of resonant planetary 
systems is very unlikely, since each resonance requires certain ratio of the semi¬ 
major axes. The more favorable scenario is the type II migration of giant planets 
being embedded in the still gas rich protoplanetary disc. Type II migration appears 
when a massive giant planet carves a gap in the gaseous protoplanetary disc prac¬ 
tically inhibiting the gas flow through the gap. In that case the planet’s semi-major 
axis is changing according to the viscous evolution of the protoplanetary disc, see 
more about the topic in (li3l). 

If the migration of two giant planets is convergent (e.g. the difference between 
their semi-major axes is decreasing) the phenomenon of the resonant capture will 
occur between them, and the two planets can migrate very close to their host star. 
The efficiency of migration is excellently demonstrated by hydrodynamic simula¬ 
tions modelling the formation of the resonant system around the star GJ 876 ll23l . 
There are other planetary systems in which the giant planets reached the 2:1 MMR 
through type II migration such as HD 128311 f38l . and HD 73526 j39l . 

The detection of giant planets is based on the radial velocity method. To calculate 
the orbital elements of the planets of a multiple system is not an easy task. Although 
there are well-known and widely used algorithms to provide reliable orbital fits, the 
orbital elements obtained not always result in a stable configuration for the planetary 
system. Regarding the resonant system of giant planets around HD 73526, ED Pub¬ 
lished orbital elements and planetary/stellar masses which resulted in stable orbits 
of the giant planets over 1 million years. On the other hand, the orbits are chaotic, 
as was clearly visible from numerical integrations of the three-body problem using 
the given elements as initial conditions. Since chaotic behavior may be uncommon 











Fig. 18 Stability maps calculated around one set of orbital elements given in |39j It can be seen 
from the stability maps that the orbital elements (marked by +) are embedded well in the stable 
region marked by white colour. 


among the resonant extrasolar planetary systems, and may not guarantee the stabil¬ 
ity of the giant planets for the whole lifetime of the system (being certainly longer 
than 1 million years), we searched for regular orbital solutions for the giant planets 
as well j39l . As a first attempt to study the degree of the chaoticity, we mapped the 
parameter space around the solution of E2- We have calculated the stability prop¬ 
erties of the a\ — ci 2 , e\ — e 2 , M\ — M 2 , and (±j\ — Uh_ parameter planes, where a is the 
semi-major axis, e is the eccentricity, M is the mean anomaly and <77 is the longitude 
of periastron of one of the giant planets. 

In Fig. 17 the stability stmctures of the parameter planes for the semi-major axis 
and the eccentricities are displayed. During the calculation of a particular parameter 
plane the other orbital data have been kept fixed to their original values. On each 
parameter plane the stable regions are displayed by white, the weakly chaotic re¬ 
gions by grey, and the strongly chaotic regions by black colors. The values of the 
corresponding orbital data are marked on each parameter plane. By studying the 
stability maps, one can see that the orbital elements given by B7l are located in a 
weakly chaotic region, which explains the irregular behavior of the planetary ec¬ 
centricities. We again stress that this does not automatically imply the instability of 
their fit, however by using these orbital data the system yields chaotic behavior and 
can be destabilized in longer timescales. Studying the stability maps it can also be 
conlcuded that the fit cannot be easily improved by the simple change of one of the 
orbital elements. The parameter plane is almost entirely weakly chaotic, there ex¬ 
ists only a narrow strip of ordered behavior. After obtaining completely new sets of 
orbital elements using the Systemic Console ( 11301 ) the stability maps around these 
fits were recalculated. In Fig. p~8] parameter planes for the semi-major axis and the 
eccentricities are displayed, now around one of the stable orbital solutions. It can be 
clearly seen that the new orbital solution is well inside the region for ordered motion, 
which provides stability for the whole lifetime of the system. Based on the above 
example, it can be concluded that the RLI performs well in stability investigations 
of extrasolar planetary systems. 
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Fig. 19 The 8:5 inner MMR (color online). Left panel: the two regions of libration as displayed by 
the libration amplitude of the resonant variable. The exact location of the MMR is at a r = 0.7310 
in normalized units. Right panel: the RLI map around the resonance, in which the two regions 
of libration can be clearly seen. The continuous curve denotes the unit apocentre distance, while 
the dotted and dashed curves the unit apocentre distance decreased by 1, and 1.5 Hills radius, 
respectively. 


3.3 Application of the RLI to study libration inside high order 
MMRs 


The most recent application of the RLI in dynamical astronomy has been presented 
by 02] investigating the dynamical structure of high order resonances in the elliptic 
restricted three-body problem. This study has relevance in dynamics of Kuiper-belt 
objects. Moreover, it proved through numerical integration of a large set of orbits 
that the RLI excellently indicates the libration of the resonant angle inside a MMR. 
In what follows, first we present the results of m obtained for the 3rd order 8 : 5 
inner resonance, in which case according to Eq. 0 the resonant variable is 

Q = 8A' — 5A — 3GJ, (7) 


where A and GJ are the mean orbital longitude, and longitude of the pericenter of the 
inner body, while A' is the mean orbital longitude of the outer body. 

A simple, but computationally demanding way to map the neighbourhood a 
MMR on the a — e plane is to calculate the libration amplitude of 0 for orbits whose 
initial semi-major axis and eccentricity values are taken from a grid, and the test 
particle has been started from pericentre. Such calculation can be seen on the left 
panel of Fig. 19 Using the same grid resolution and initial conditions the RLI val¬ 
ues were also calculated, see the right panel of Fig. 19 The exact location of the 
resonance (when Q = 0°) is marked with the vertical dashed line in the left panel of 
Fig. [19] The colour code corresponds to the value of the libration’s amplitude: the 
darker the colour, the smaller the amplitude. 

The 8 : 5 MMR has an interesting structure, there is a libration for low values of 
the eccentricities e < 0.2, and also for very high values 0.75 < e < 0.85. Although 
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Fig. 20 The 8:5 outer MMR (color online). Left panel: the two regions of libration displayed by 
the libration amplitude of the resonant variable. The exact location of the MMR is at a r = 1.3680 in 
normalized units. Right panel: the RLI map around the resonance, in which the regions of libration 
are clearly visible. The continuous curve indicates the unit pericentre distance, while the dotted 
curve the unit pericentre distance increased by one Hills radius. The dash-dotted curve shows the 
place of Uranus-crossing. The filled circles mark circulating TNOs. 


the picture obtained by the RLI is more detailed, also showing some other neighbour 
MMRs, gives back the region of libration very well. 

The RLI has also been applied to study the high order outer MMRs, in which case 
the massless body’s orbit is outside the massive planet’s orbit. In the study of m 
the outer resonances of the Sun-Neptune system have been studied in two different 
models. These are the restricted three-body problem, in which only the gravitational 
effects of the Sun and Neptune were included, and also in a model in which the four 
giant planets were also taken into account. We will present here the results obtained 
for two different MMRs, namely the 8:5 and 7:3 outer resonances in the model of 
the restricted three-body problem. We note that the resonant variable of an outer 
(p + q) : p MMR is given by Eq. ||6}. 

Similarly to the 8 : 5 inner MMR, the neighbourhood of the exact 8 : 5 outer 
MMR ( a r = 1.3680 in normalized units) has been mapped on a dense grid of the a — 
e plane (the test particle has been started from apocentre) and the libration regions 
are marked with different shades corresponding to the amplitude of 9. There are 
two regions of libration in this resonance located at relatively high eccentricities 
(0.35 < e < 0.7, 0.75 < e ), see the left panel of Fig. |20| Using the same a — e 
grid the RLI values have also been calculated, see the right panel of Fig. 20 At 


the location of the exact resonance there is a nearly vertical dark strip indicating 
weak chaotic behaviour, and also the fact that in this configuration the 8 : 5 MMR 
is not protective for low values of the test particle’s eccentricity. The origin of the 
strong chaotic behaviour is due to the crossing of the orbits of the test particle and 
the perturbing body (Neptune in this case). There are different lines plotted to the 
RLI map. The continuous curve is the pericentre distance, the dotted curve is the 
pericentre distance increased by the Hill’s radius of Neptune, and finally, the dash- 
dotted curve is the place of Uranus crossing. This latter curve may indicate that 
the high eccentricity libration regions of Fig. [20] being present in the model of the 
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Fig. 21 The 7:4 MMR (color online). Left panel: the two regions of libration displayed by the 
libration of the resonant variable. The exact location of the resonance is at a r = 1.4522. Right 
panel: the RLI map around the resonance. The triangles mark librating, while the filled circles 
circulating TNOs. For the description of the curves see the previous figures. 


restricted three-body problem, might vanish in more advanced models including the 
planet Uranus, for instance. Finally, we roughly compare the 8 : 5 outer resonance 
with the behaviour of some of the existing TNOs having inclination i < 25°, see the 
black filled circles for the corresponding a and e values of these objects. All of the 
TNOs displayed are circulating, and really they occupy the low eccentricity regions 
of the a — e plane. We note, however, that the a and e values of the TNOs have 
different epochs, so their positions does not reflect the actual state of the system. On 
the other hand, it can be clearly seen in Fig.[20]that the TNOs at the 8 : 5 MMR have 
circulating resonance variable. 

In order to have a more complete picture of the high order outer resonances, we 
also summarize the case when a MMR has a protective character, and the resonant 
variable of bodies lying in its vicinity can both librate and circulate. A good can¬ 
didate for this purpose is the 7 : 4 3rd order outer MMR. The exact resonance is at 
a r = 1.4522 (normalized units). Similarly to the 8 : 5 outer MMR, this resonance 
also has two regions of libration, but in this case the lower region of libration al¬ 
lows libration of test particles having low eccentricities, see Fig. 21 The right panel 
shows the dynamical structure of the resonance, and also TNOs found at this reso¬ 
nance. Most of them have circulating 0, but there are a few of them in the librating 
region, too. Studying the right panel of Fig. 21 one can see that the librating TNOs 
are clearly below the Neptune-crossing line, while the region of libration extends 
above it. 

As an overall conclusion we can state that the RLI is a very reliable tool in de¬ 
tecting the positions of the lower and higher order MMRs being present in various 
planetary systems. 
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4 Discussion and summary 

In this chapter we summarize the basic properties of the recently introduced chaos 
detection method, the RLI. The RLI is based on the time evolution of the infinitesi¬ 
mally small tangent vector to the orbit, which is provided by solving numerically the 
variational equations. Hence, the RLI belongs to the family of the so-called varia¬ 
tional indicators. Although the definition of the RLI is based on that of the LI, in this 
review we give evidence that the distinction between regular and chaotic motion is 
much clearer with the RLI, which makes it a more reliable alternative than the LI. 

According to the comparative study with some wide-spread variational indi¬ 
cators, the RLI shows convincing performances in the experiments and consider¬ 
ably improves the performances of the classical LI. In generality, indicators like the 
FLI/OFLI or the MEGNO (actually there is a strong relationship between both in¬ 
dicators, see ED for further details) are usually believed to be better options for a 
general analysis of the structure of the phase space. Therefore, our study reinforces 
the fact that the RLI can also be used as an alternative technique, which operates 
with reasonable computing times to make conclusive pictures of the dynamics de¬ 
spite the complexity of the problem. 

Based on the comparative work presented in Section [2j we can summarize both 
the advantages and disadvantages of the RLI. In what follows we list its favorable 
properties/advantages, and also add that in which dynamical systems have been done 
the corresponding simulations. 

• Having determined a reliable threshold in the Henon-Heiles model, the RLI (and 
the FLI/OFLI) shows the best approximation rates in the ordered/chaotic regions. 
We note that the methods SALI and GALI also estimate the true percentage of 
the ordered/chaotic orbits but with a slight slower way. 

• Comparing to the other chaos indicators also in the Henon-Heiles model, the RLI 
detects much faster the orbits from the large chaotic sea (e.g. the “c-cs” orbit), 
than the other indicators. 

• Studying a 4D symplectic mapping the RLI have been compared to the indi¬ 
cator MEGNO(2,0). In this case both the RLI and the MEGNO(2,0) reveal the 
fine structure of the phase space very accurately (much better than the LI, for 
instance). As a result of this experiment we also conclude that the RLI is a very 
effective tool in the characterization of a large array of initial conditions. 

• In a complex 3D potential resembling a triaxial galactic halo (the so called NFW 
model), the RLI (together with the MEGNO, the OFLI, and the GALI 5 ) identifies 
the chaotic orbits within a Hubble time (~ 13 Gyrs). These indicators show that 
chaotic orbits can be identified within a physically meaningful time (i.e. the age 
of the Universe), which is important when studying the dynamics of a galaxy. 

On the other hand, when applying the RLI one should be aware that: 

• Among the presented CIs, in the Henon-Heiles model the RLI identifies the so 
called “c-sl” orbit in the slowest way. 
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• In the study performed in the 4D symplectic mapping, the RLI cannot really 
distinguish between chaotic and sticky orbits. This is a disadvantage if one is 
interested in detecting the sticky orbits. On the other hand, if we are interested 
in detecting all chaotic orbits (including sticky orbits, as well), the application of 
the RLI might be useful. 

• Regarding the time needed to calculate the RLI, we can conclude that it is not 
the least time consuming indicator. The FLI/OFLI, the MEGNO, the SALI and 
the SD might be more desirable options if the computation time really matters. 
We note, however, that with the current generation of fast computers this option 
became less important. 

In the last section of the current work we summarize the application of the RLI to 
planetary systems, which is its major application area. These studies include the 
development of a stability catalogue of hypothetical terrestrial exoplanets in ex¬ 
trasolar planetary systems, stability studies of resonant planetary systems and the 
investigation of high order mean motion resonances having relevance in studying 
the dynamics of the Kuiper-belt objects. We find that the RLI is an efficient and 
reliable numerical tool to map and characterize the dynamical structure of various 
mean motion resonances. 

We note that the RLI (together with the SALI) has also been applied to map the 
stability regions of the Caledonian symmetric four-body problem, ( ll48ll ). Since the 
preceeding studies are mainly related to detecting the chaotic behaviour occured 
near to resonances, the study on the Caledonian restricted four-body problem does 
not really belong to this line, thus to keep the lenght of the present study tractable, 
we omited its presentation here. 

We would also like to remark that the very simple computation of the RLI from 
the widespread well-known LI and its better performances reported in very differ¬ 
ent scenarios, make the RLI a serious candidate to replace the LI in a variety of 
fields, and not only in dynamical astronomy. For instance, in a paper published in 
a journal of Chemical Physics, the RLI is used as the default chaos indicator in the 
Lyapunov weighted path ensemble method. One of the capabilities of the method 
is to identify pathways connecting stable states which are relevant in the context of 
activated chemical reactions (see im 

Finally, the reliability of the RLI as a chaos indicator has been strongly demon¬ 
strated throughout this study, and as a result, the choice of the RLI to analyze a 
general dynamical system is well-founded. 
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